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THERMOACOUSTIC TOMOGRAPHY WITH AN ARBITRARY ELLIPTIC OPERATOR 

MICHAEL V. KLIBANOV* 

Abstract. Thermoacoustic tomography is a term for the inverse problem of determining of one of initial conditions of a 
hyperbolic equation from boundary measurements. In the past publications both stability estimates and convergent numerical 
methods for this problem were obtained only under some restrictive conditions imposed on the principal part of the elliptic 
operator. In this paper logarithmic stability estimates are obatined for an arbitrary variable principal part of that operator. 
Convergence of the Quasi-Reversibility Method to the exact solution is also established for this case. Both complete and 
incomplete data collection cases arc considered. 
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\Q ■ 1. Introduction. The goal of this paper is to show that logarithmic stability estimates as well as 

convergent numerical methods for the inverse problem of determining an initial condition in a general hy- 
perbolic PDE of the second order can be obtained without any restrictions on its coefficients, except of 
<— j ] some natural ones. In all previous publications on this topic the principal part of the elliptic operator was 

Oh subjected to some restrictive conditions. Naturally, our stability estimates imply uniqueness. Both complete 

and incomplete data collection cases are considered. We assume here that the data are given on the infinite 
time interval t £ (0, oo) . Second and third Remarks 2.1 (section 2) justify this assumption. For brevity, we 
leave for possible future publications the finest assumptions, like, e.g. the minimal smoothness, etc. 

In thermoacoustic tomography (TAT) a short radio frequency pulse is sent in a biological tissue [2[ [TO]. 
Some energy is absorbed. It is well known that malignant legions absorb more energy than healthy ones. 
Then the tissue expands and radiates a pressure wave which is the solution of the following Cauchy problem 
>' 

ii ii = c 2 (i)Ati,ieM 3 ,i>0, (1.1) 
u (x, 0) = / (x) , u t (x, 0) = 0. (1.2) 



The function u (x, t) is measured by transducers at certain locations either at the boundary of the medium of 
interest or outside of this medium. The function / (x) characterizes the absorption of the medium. Hence, if 
one would know the function / (x), then one would know locations of malignant spots. The inverse problem 
consists in determining / (x) using those measurements. 

Both stability estimates and convergent numerical methods for the problem of determining the initial 
condition / in (jl.2| are currently known only under some restrictive conditions imposed on the coefficient 
c (x) (subsection 1.2). In addition, except of the case c(x) = 1 in [21] . those numerical methods are known 
5_i " only for the case of complete data collection, i.e. when boundary measurements are given at the entire 

boundary of the domain of interest. 

First, we apply a well known analog of the Laplace transform to obtain a similar inverse problem for 
a parabolic PDE. Next, previous results of the author [HI [17] are used. In the complete data case the 
logarithmic stability estimate follows from [16] . In the case when the data are given on a hyperplane, 
we significantly modify the proof of Theorem 1 of |17j . More precisely, we prove our logarithmic stability 
estimate for an integral inequality rather than for the parabolic PDE. We need this generalization to establish 
convergence rate of our numerical method. Results of both publications [16[ I17j were obtained via Carleman 
estimates. In particular, a quite technical non-standard Carleman estimate was derived in [17] . see Lemma 
2.1 in section 2. We refer to [31] for another logarithmic stability estimate of the initial condition of a 
parabolic equation with the self-adjoint operator L in a finite domain. A Carleman estimate was also used 
in this reference. An interesting feature of [31] is that observations are performed on an internal subdomain 
for times t € (r, T) where r > 0. In addition, a numerical method was developed in [31] . 
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1.1. Statements of inverse problems. Let SI C {x\ > 0} be a bounded domain with the boundary 
dn G C 3 . Let T > 0. Denote 

Q T = O x (0, T) , S T = dfl x (0, T) , P = {xi = 0} , P T = P x (0, T) . 

Let > be an integer and a G (0, 1). Below C k+a , C 2k+a ' k+a l 2 are Holder spaces. Consider the elliptic 
operator L of the second order with its principal part Lq, 

n n 

Lu = dij (x) u XiX] + ^2 b j (%) u X] + b (x) u, x G R", (1.3) 

i,j=l 3=1 
n 

L u = ^ a i-j (x)u XiXj , (1.4) 
i,j=i 

a hJ £ C k+a (W 1 ) n C 1 (R n ) , bj,bo € C k+a (R n ) , k > 2, a € (0, 1) , (1.5) 

71 

Mi M 2 < X! a« j (z)^ ? ?j ■< M2 l^l 2 >Vx,7y S M"; /x l5 /x 2 = const. > 0. (1.6) 

Let the function / (x) be such that 

f eC p (M. n ),p>3J(x) =0,xGM"\O. (1.7) 
Consider the following Cauchy problem 

Uu = Lu, x G M™, t G (0, co) , (1.8) 
u(x,0) = /(x),tt t (x,0) = 0. (1.9) 

We use everywhere below the following assumption. 

Assumption. We assume that integers k > 2,p > 3 in (|1.5p . (|1.7p . coefficients of the operator L and 
the initial condition / are such that there exists unique solution u G C 3 (R™ x [0, T]) , VT > of the problem 
(|T8j) . (fl~9l> satisfying 

ll«llo«(R»x[o,ri) <Be dT ,VT>0, (1.10) 

where the constants B = B (L) > 0,d = d(L) > depend only from the coefficients of the operator L and 
an upper estimate B of the norm < B. 

Note that (|1.7p in combination with the finite speed of propagation of the solution of problem (|1.8I) . 
(UU) guarantee that the function u (x,t) has a finite support * (T) C R",Vi G (0,T) ,VT > [25]. Hence, 



C (R" x [0, T}) above is actually the space C (T) x [0, T]j . Using the classical tool of energy estimates 

[25] . one can easily find non-restrictive sufficient conditions imposed on both coefficients of the operator L 

and the function / guaranteeing the smoothness u G C 3 (R™ x [0,T]) ,VT > as well as (|1.10[) . We are not 

doing this here for brevity. We consider the following two Inverse Problems. 

Inverse Problem 1 (IP1). Suppose that conditions S 1.3)) - (T7)\ ) and Assumption hold. Let 

ueC 3 (R n x [0,T]) , VT > be the solution of the problem fO)) . IO|) . A ssume that the function f (x) 

is unknown. Determine this function, assuming that the following function tp 1 (x, t) is known 

u\ Soe =<Pi(x,t). (1.11) 

Inverse Problem 2 (IP2). Suppose that conditions $1.3\ )- (T7)\ ) are Assumption hold. Let u G 
C 3 (R™ x [0,T]),VT > be the solution of the problem 111-8]) . il.Sfy . Assume that the function f (x) is 
unknown. Determine this function, assuming that the following function ip 2 (x, t) is known 

W \x€Pco= ¥2 (pit) ■ ( L12 ) 

IP1 has complete data collection, since the function ip 1 is known at the entire boundary of the domain of 
interest ft. On the other hand, IP2 represents a special case of incomplete data collection, since ft C {x\ > 0} . 
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1.2. Brief overview of published results. TAT has attracted a significant interest in the past several 
years. We now provide a brief overview of published mathematical results for TAT. We refer to |20j for a 
review paper. Stability estimates and convergent numerical methods for an arbitrary time independent 
principal part Lq in (jl.4[) were not obtained in the past. Explicit formulas for the reconstruction of the 
function / (x) for IP1 in the case when in c = 1 arc given in a number of publications, see, e.g. 

[8j [9] [TOl [20] [21] . These formulas lead to some stability estimates as well as to numerical methods with good 
performances. 

Another approach to IP1, IP2 is via analyzing the case when both Dirichlet and Neumann data are 
given at St for IP1 and at Pt for IP2. An elementary, well known and stable procedure of deriving the 
Neumann condition from the given Dirichlet condition for both IP1 and IP2 is described in subsection 2.1 
for the parabolic PDE. A very similar procedure takes place in the hyperbolic case. Consider now IP1. Since 
a certain norm of the Neumann boundary condition at St can be estimated from the above by another 
norm of the data ip 1 (x,t) for (x,t) G St, then the problem of determining the initial condition f (x) can 
be reformulated in a slightly more general setting as the Cauchy problem for equation (|1.8[) with the lateral 
Dirichlet and Neumann data at St- This problem consists in determining the function u(x,t) inside of the 
time cylinder Qt- 

We now comment on the Lipschitz stability estimate for that Cauchy problem with lateral data for 
the particular case when initial conditions are as in (|1.9[) . Consider the even extension of the function 
u(x, t) with respect to t and do not change notations for brevity, u (x, — t) := u(x,t),t € (0,T). Let 



Qt 



fl x (-T, T),St — dflx (-T, T) . Obviously 



2 II" \s T \\m(s T ) and 



d v u 



2 \\d v u \s t \\l- 2 {St) ' wnere ®v means the normal derivative. The Lipschitz stability estimate for the Cauchy 
problem with the lateral data is 



Nli3-i(Q±) < C \\\u \ St \\ h i (St) + \\d„u \s t \\ L2{St ) 



(1.13) 



with a certain constant C > independent on the function u. Hence, the trace theorem implies the Lipschitz 
stability estimate for the function / with a different constant C, 

H"(^0)llL 2 (rr> = II/IIl 2 (J2) < c \\ u \st\\hhs t ) + R u \s t \\l 3 (s t ) 



Estimate (|1.13|) plays a fundamental role in the control theory, since it is used for proofs of exact 
controllability theorems. For the first time estimate (|1.13|) was proved in 1986 in (32] for equation ([l.ljl 
with c = 1 with the aim of applying to the control theory. However, the method of multipliers, which was 
proposed in [32] , cannot handle neither variable lower order terms of the operator L nor a variable coefficient 
c (x) . On the other hand, Carleman estimates are not sensitive to lower order terms of PDE operators and 
also can handle the case of a variable coefficient c (x) . 

In [13] the Carleman estimate was applied for the first time to obtain (|1.13j) . In [13] (|1.13j) was proved 
for the case of the hyperbolic equation (|1.8|) with L = A+ (variable lower order terms). Next, the result of 
[T5] was extended in [T^] [T3] to a more general case of the hyperbolic inequality 



\u tt ~Au\ <A[|V«| 



I/O in Qt, 



(1.14) 



where A = const. > and / G (Qt) ■ Although in publications JT2l [13] [15] c = 1, it is clear from them 
that the key idea is in applying the Carleman estimate, while a specific form of the principal part of the 
hyperbolic operator is less important. This thought is reflected in the proof of Theorem 3.4.8 of the book 
|11) . Thus, the Lipschitz stability estimate (|1.13l) for the variable coefficient c(x) was obtained in section 
2.4 of the book [18] as well as in [6]. In particular, in [18] the case of an analog of the hyperbolic inequality 
(|1.14j) was considered, where \u u — Au| was replaced with |c -2 (x) u tt — Au\. The idea of [12] was used in 
the control theory in, e.g. |27 [ I28 j . 

In the case of parabolic and elliptic operators, Carleman estimates are known for rather arbitrary variable 
principal parts [TT] [TH] [3D]. On the other hand, it is well known that in the hyperbolic case the Carleman 
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estimate can be effectively analytically verified for a generic operator d 2 — Lq only if L = c 2 (x) A, and a 
condition like 



(a; -so, V (c~ 2 (x))) > 0, Vx e SI 



(1.15) 



holds. In (|1.15p Xo is a certain point and (, ) is the scalar product in M™. This is the reason why the 
above mentioned Lipschitz stability estimates were established only using assumptions like the one in (|1.15l) . 
Clearly, (|1.15p holds for c = const. ^ 0. See, e.g. Theorem 1.10.2 in [5] for the proof of the Carleman 
estimate with condition (|1.15[) . A more general case of condition (|1.15[) can be found in Theorem 3.4.1 
of [11]. The second way of proving Lipschitz stability estimates is via imposing some conditions of the 
Riemannian geometry on coefficients of the operator L [H |26l [34j [35] [36] . Publications [26[ EH ES] use 
some modifications of the idea of |12 [ 113 ) : in particular, the case of hyperbolic inequality was considered in 
[34] . Unlike (|1.15|) . conditions of the Riemannian geometry cannot be effectively analytically verified for an 
operator Lq with generic coefficients, e.g. Lq ~ c 2 (x) A. A slight variation of (|1.15[) guarantees the non- 
trapping condition, see formula (3.24) in [53]. Uniqueness theorems for TAT were also obtained in [2"1 HUl 156] 
for the case ([L~T]) . ([L2]) . 

In addition, to the Lipschitz stability, the Quasi-Reversibility Method (QRM) for the above mentioned 
Cauchy problem with the lateral data was developed in [T5J and numerically tested in JS] Q3J Q]5] . We refer 
to [22] for the originating work on QRM. The convergence of the QRM solution to the exact solution was 
proven on the basis of the above Lipschitz stability results. Numerical testing has consistently demonstrated 
a high degree of robustness. In particular, accurate results were obtained in [19j with up to 50% noise in 
the data. Some other numerical methods were proposed in [2] [36] . Convergence of all numerical methods 
mentioned in this paragraph was proven only for the complete data collection case of IP1 with Lq = c 2 (x) A 
and under some restrictive conditions imposed on the function c (x). 

2. Logarithmic Stability. 

2.1. Transformation. First, we consider the following well known Laplace-like transformation [T51I50"] . 
which transforms the hyperbolic Cauchy problem in a similar parabolic Cauchy problem, 



Cg = g(t) = 



1 



cxp 



4/ 



g (r) dr. 



(2.1) 



The transformation (|2.1[) is an analog of the Laplace transform, and it is one-to-one. It is valid for, e.g. all 
functions g e C[0,oo) which satisfy \g{t)\ < A g e kat , where A g and k g are positive constants depending on 
g. It follows from (|1.10[) that the solution u (x, t) of the problem (|1.8[) . (|1 .9[) satisfies this condition together 
with its derivatives up to the third order. Obviously 



Hence, 



d_ 

at 



i 



exp 



T 



dr 2 



1 



exp 



T 

4f 



C (g") = g' (t) , V.9 e C 2 [0, oo) such that g' (0) = 0. 



(2.2) 



Changing variables in f|2 . 1 [) r <^=> z,r/2^/t := z, we obtain lim t ^ + .9 (t) =5(0)- Denote v := Cu. It follows 
from (fTTTO]) and ([22]) that 

v G C 2+Q ' 1+a/2 (R n x [0, T] ) , Va G (0, 1) , VT > 0. (2.3) 
By (|1.8p . (|1.9p and (|2.3p the function v (x,t) is the solution of the following parabolic Cauchy problem 



v t = Lv,x G W\t> 0, 
u(x,0) = /(x). 
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(2.4) 
(2.5) 



We refer here to the well known uniqueness result for the solution v £ (j2+a,i+a/2 x [0,T]) , VT > of 
the problem ([£3]) . (j2~5l> [SJ. 

Below we work only with the function As to this function, we set everywhere below T := 1 for the 
sake of definiteness. Denote 

£</>i = (2;, t) =v\ Sl , Cip 2 = Tp 2 (x, t) = v | Pl . (2.6) 

Then 

- e c 2 +Q ,l+«/2 fij ^ g C ,2 +a ,l+a/2 (p^ _ (2 7) 

Let 

"01 (^*) = 9»« lsi,V , 2 OM) = 9 xi v \ Pl . (2.8) 

By Theorem 5.2 of Chapter IV of M, (Ti~TUj) and (JUJl-dl^) there exist numbers C (Q, L) ,C (P, L) > 
depending only on listed parameters such that 

||^l|| c i+ Q ,o/2(5 1 ) < C (n,L) H^ill^+c.i+a/a^) , (2.9) 

||V , 2|| C l + c,, Q /2 ( --p 1 ) < C{P,L) 1 II c2 + Q ,l + cl/ 2(p 1 ) ■ (2.10) 

We now describe an elementary and well known procedure of finding the normal derivative of the function 
v either at Si (in the case of IP1) or at Pi (in the case of IP2). In the case of IP1 we solve the initial 
boundary value problem for equation (|2.4|) for (x,t) G (K"\f2) x (0,1) with the zero initial condition in 
IR n \S! (because of (jl -7|) ) and the Dirichlet boundary condition v \s! = ~<Pi- Then we uniquely find the normal 
derivative d u v \s 1 = Similarly, in the case of IP2, we uniquely find the Neumann boundary condition 
9 X i v \pi = Vv Estimates (|2.9p . (|2 . 10[) ensure the stability of this procedure. 

Therefore, each problem IP1, IP2 is replaced with a problem for the parabolic PDE (|2.4j) with the lateral 
Cauchy data. These data are given at Si for IP1 and at Pi for IP2. Uniqueness of the solution of each of 
these parabolic inverse problems follows from standard theorems about uniqueness of the continuation of 
solutions of parabolic PDEs with the data at the lateral surface |11[ [151 [50] . 

In stability estimates one is usually interested to see how the solution varies for a small variation of the 
input data. Therefore, following (|1.10[) . and (|1.12[) . we assume that in the case of IP1 

IMIc3(S T ) <<5e rfT ,VT>0, (2.11) 

and in the case of IP2 

\\<P2\\c(p T ) <5e dT ,VT>0, (2.12) 

where 5 £ (0, 1) is a sufficiently small number. Note that it is not necessary that S = B, where B is the 
number from (|1.10p . Indeed, while the number B in (|1.10[) is not assumed to be sufficiently small and is 
involved in the estimate of the norm IM| C 3( K „ X [ t]) > in the entire space, the number 5 is a part of 
the estimate of the norm of the boundary data for either of above inverse problems. Using (|2.1[) . (|2.2I) and 
(|2l ^ - ([2TT2|) . we obtain 

||^i|| C 2+ 0! ,i+o,/2(g 1 ) + ||'0i|| c i+ Q , Q /2(5 1 ) < C (0, L, d) S, (2.13) 

H^llcs+o.i+a/^F!) + ||V'2|lci+<w(F 1 ) <C(P,L,d)8, (2.14) 

where constants C (fi, L,d) , C (P, L,d) > depend only on listed parameters. It follows from (|2.13[) that 
with a different constant C := C (fl, L, d) > 



(2.15) 



Remarks 2.1. 

1. The number 6 can be viewed as an upper estimate of the level of error in the data tp 1 ,tp 2 . Hence, 
Theorems 2.1, 2.2 below address the question of estimating variations of the solution / of either IP1 or IP2 
via the upper estimate of the level of error in the data. 

2. Since the kernel of the transform C decays rapidly with r — > oo, then the condition t £ (0,oo) in 
(jl.lip . (|1.12p is not a serious restriction from the applied standpoint. In addition, if having the data in 
(jl.ll|) . (|1.12j) only on a finite time interval t € (0,T) and knowing an upper estimate of a norm of the 
function / in (|1.9p . one can estimate the error in the integral (|2.ip when integrating over r <G (T, oo) . This 
error will be small if either T is large or t is small in (|2.ip . (|2.4p . Next, this error can be incorporated in the 
stability estimates of theorems of this section. 

3. Another argument about t £ (0, oo) comes from the recent experience of the author of working 
with time resolved real data for wave processes [5]. The author has learned that almost all time resolved 
experimental data for wave processes in non-attenuating media are highly oscillatory due to some unknown 
processes in measurement devices, see graphs of those data in these references. Because of high oscillations, 
these data are not governed by a hyperbolic PDE even for the case of the free space, where the wave equation 
is supposed to work (see the graphs of experimental data in chapters 5 and 6 of [5]). Therefore, the first step 
to make the inverse algorithm work was to preprocess the experimental data via a new data preprocessing 
procedure. This procedure uses only a small portion of the real data and immerses it in a specially processed 
data for the uniform medium. Since the case of the uniform medium can be solved analytically, then there 
is no problem to know the immersed data for all t £ (0,oo) . Since accurate imaging results were obtained 
in [5] for the case of blind experimental data, then that data preprocessing procedure was unbiased. 

2.2. Logarithmic stability estimate for Inverse Problem 1. To prove convergence of the QRM 
(Theorem 3.1), it is convenient to consider a parabolic inequality in the integral form, which is more general 
than equation (|2.4p . Consider the function w £ C 2,1 (Qi) satisfying the following inequality 

J (w t - Lwf dxdt < K 2 ,K = const. > 0. (2.16) 

Qi 

Theorem 2.1. Let conditions U.3\) - [J~S\) be fulfilled. Let the function w £ C 2,1 (Qi) satisfies inequality 
V2~W\) . Denote 

g(x) =w(x,0),p (x,t) =w |si,£i (x,t) =d v w 

F = m\HH Sl ) + m\L 2{Sl )+ K - (2-i7) 

Assume that an upper bound C\ = const. > for the norm ||V<7|| ia (n) is known, 

lL 2 (H) — II v .yilL 2 (fi) 



lis* 



"Vff|| £ilfm < C x . (2.18) 



v i=l 



Then there exist a constant M = M (L, Q) > and a sufficiently small number Sq = 8q (L, f2, Cx) £ (0, 1) , 
both dependent only on listed parameters, such that if F £ (0,<5o); then the following logarithmic stability 
estimate is valid 

\\g\\j lm < MCl . (2.19) 

Ln particular, in the case of IP1, let Assumption holds and \1. 7| ), V2.ll]) be valid. Suppose that the number 
5 in &2.11]) is so small that CS £ (0, do), where C = C (Sl,L,d) > is the number in i2.15)) . Also, assume 
that the upper bound C\ of the norm ||V/|| i2 /™ is given, 

l|V/|| ia(n) < Cx. (2.20) 
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Then 




< 



MCi 



(2-21) 



In (CS) 



1 



Proof. In this proof M ~ M (L, fi) > denotes a generic positive constant depending only on L, il. 
First, we prove (|2.19|) . Let x £ (0, 1) be an arbitrary number. Then it follows from Theorem 2 of [16] that 
there exists a constant r — r (L, il, x) £ (0, 1) such that 



as long as F € (0,1). We can fix x via, e.g. setting x := 1/2. It is clear therefore that there exists a 
sufficiently small number S = S (L, Q, C\) £ (0, 1) such that if F £ (0, S ) , then ([2~22]) implies ([2TT9]) . 

We now prove (|2.21[) . It follows from (I2.4[) that (|2.16|) holds for the function w :— v with K = 0. As 
it was shown above, (|2.15[) follows from (|2.11l) . Hence, using (|2.15[) and (|2.17[) . we obtain F < CS. Hence, 
([2~2l"jl follows from (j2~19| . □ 

Remark 2.2. Estimates (|2.19|) . (|2.21[) are the so-called "conditional stability estimates", which is often 
the case in ill-posed problems [5j [30] [37] . For another example we refer to Holder stability estimates for 
solutions of some ill-posed problems for PDEs, see, e.g. [11] [T8] [30]. The knowledge of the upper bound 
C\ for the gradient in (|2.18j) . ()2.20j) corresponds well with the Tikhonov concept of compact sets as sets 
of "admissible" solutions of ill-posed problems [3] [5] [7] [30] [37]. Indeed, since by (|1.7[) / |an= 0, then 
||/||^, 2 (f2) < R ||V/|| L2 (f 2 -) < RCi, where the constant R > depends only on the domain H. Thus, in this 
case the function / belongs to a compact set in L2 (fi) , and this set is determined by the constant C±. 

2.3. Logarithmic stability estimate for Inverse Problem 2. The logarithmic stability estimate 
of the paper |17j in the infinite domain was obtained for the case of the pointwise inequality 



where the operator Lq is defined in (|1.4[) . However, to prove convergence of the numerical method of section 
3, we need to estimate the initial condition for the case of the integral inequality, like the one in (|2. 16[1 . The 
Carleman estimate of |17| is not a standard one. Indeed, unlike the standard Carleman estimate for the 
parabolic operator [21 UBl ED] , the integration domain of fTT] is a part of the strip {|f — e| < re, r £ (0, 1)} , 
and that Carleman estimate does not break when e-)0 + . 

There are two main differences between Theorem 2.2 (below) and Theorem 1 of |17| . First, we work 
now with the integral inequality instead of the pointwise inequality (|2.23[) of |17j . Second, it is assumed in 
[T7] that the inequality (|2.23p is valid in 9 x (0, T) , where C R™ is an unbounded domain. It is also 
assumed that the Dirichlet boundary condition v laex(o.T)— 0. Unlike this, Theorem 2.2 does not use the 
assumption about the knowledge of this Dirichlet boundary condition. 

Denote x = (X2, x n ) . Below we specify numbers 1/4,1/2,3/4 for brevity only. In fact, some other 
numbers, respectively rj 1 < rj 2 < ?y 3 < 1 from the interval (0, 1) can be used. Changing variables (x',t') = 
{■sfcx, dt) with an appropriate constant c > and keeping the same notations for new variables for brevity, 
we obtain that 




(2.22) 



\v t - L v\ < A (|Vu| + M) , A = const. > 0, 



(2.23) 




(2.24) 



Let e £ (0, 1) be a sufficiently small number 



Consider the following functions ip (x, t) , ip (x, t) , 




(2.25) 




(2.26) 
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where v > 1 is a large parameter which will be defined later. The function ip (x, t) is the Carleman Weight 
Function (CWF) in the Carleman estimate of Lemma 3.1. The main difference between <p (x, t) in (|2.26|) and 
the standard CWF for the parabolic operator [SJ [T5J [3D] is that the small parameter e is involved in both 
functions ip(x,t) and ip(x,t). Denote 

(2.27) 
(2.28) 



G 3 / 4 = ]{x,t) : i> (x,t) < -,x x > \ , 



G1/2 = \ (x,t) :ip{x,t) < -,xi > 



Using (|2~2"2|) - ([2~2"g|) . we obtain 

G1/2 C G3/4, ip 2 (x, t) > exp 



0H-1 



in G1/2 



G3/4 C { \t-e\ < 



V2 



C{t 6(0,1)}, 



(2.29) 

(2.30) 
(2.31) 
(2.32) 

In (|2.3ip RG1/2 and RG3/4 are orthogonal projections of domain G1/2 and G3/4 respectively on the hypcr- 
planc {t = 0} . The same notation RH is kept below for the projection of any other domain H C [R" x (0, 1)] 
on the hyperplane {t = 0} . Denote 



n c RG 1/2 c i?G 3/4 , 

<9G 3 / 4 = 9iG 3/4 U d 2 G 3/4 , <9iG 3 / 4 = {xi = 0} n G 3/4 , <9 2 G 3/4 = <{ ip (x, t) = ^, xi > 



(2.33) 
(2.34) 

(2.35) 



$={(x,t): Xl G (0,1), (as, 13, -.in) G HMf'Ste (0,1)}, 
8 1 $ = inP={(z l t) : xi =0,s e (-l,l) n ~\t e (0,1)} . 

By (237) . d£30| and (l2"52l-(l2"Ml) 

diG 3/4 C 

Recall that (|2. 12[) implies (|2.14[) . Hence, assuming that (|2.12[) holds and using (|2.34[) . we derive, similarly 
with the above derivation of (|2.15j) from (|2.11j) . (|2.13|) . that there exists a constant G = C (P, L, <f>,6) > 
such that 

H^allHi( ft ») + l|Vy L2(ai<t> ) < GJ. (2.36) 

Everywhere below G = G (£0, -RG3/4) > and Mi = Mi (L, <I>) > denote different positive constants 
depending only on listed parameters. The following lemma is follows immediately from Theorem 2 of |17j 
and (2J2J. 

Lemma 2.1. Lef coefficients of the operator Lq in l{1.4\ ) satisfy conditions hl.5\) . il.6\) . Then there exist 
a sufficiently large constant v$ = Vq (Lq, -RG3/4) > 1 and a sufficiently small number Eo = Eq (Lo, RG3/4) E 
(0, 1) , both dependent only on Lq and PG3/4, such that the following Carleman estimate holds 



Cv 3 



■ exp 



2 • 4 V 



u 2 + |Vu| 2 + u 2 



dxdt 



Cv 3 fA 
' .3 



3iG 3 / 4 

~2 /4 



exp 



> G 



v 



e V3 

4 



(u 2 + |Vu| 2 + u 2 ^j da + / (lit — Lqu) 2 ip 2 (a;, t) efedi 



9 2 G 3 



|Vw| z + ~^i> u 9 (x,t)dxdt,\fu > v , Ve G (0,e o ),Vu G G 



2,1 



G 



3/4; 



? 3/4 



By (|2.30[) Lemma 2.1 provides the Carleman estimate in the narrow strip ||£ — e| < e/%/2} . At the same 
time, it is also important in numerical studies of the QRM to estimate its solution in a not narrow strip. 
This can be done via the standard Carleman estimate. Therefore, we introduce now notations, which are 
similar with (|2.25[) - (|2.32[) . except that a narrow strip with respect to t is not used. Let 



'(x,t)=x 1 + \x\'+ [*-£) + 



£(x,t) = exp(X0~ v ) , 
where A > 1 is a large parameter which is chosen later. Denote 

D3/4, = |(x,t) : 0(x t t) < > o| , 



D1/2 = \ (x,t) :9(x,t) < -,xi > 



3D 



3/4 



diD 3/4 U d 2 D 3/4 , d x D 3 ^ = {xi = 0} n D 3/4 , d 2 D 3/i = \ V> (x, t) = -, x x > 



Using P^gjl . (|OT|) . (j2~M|) . ([237)1 and (j2~3g|) - ([2~4Tj) . we obtain 

c RD 3 / 4 = RG 3 / 4l 



D3/4 c 
Ai/4 c $ 



£ - 



< - } c {t e (0, l)} , 



(2.37) 
(2.38) 

(2.39) 
(2.40) 
(2.41) 

(2.42) 
(2.43) 
(2.44) 



Lemma 2.2 follows from the Carleman estimate for the parabolic operator of Lemma 3 of §1 of Chapter 4 of 
the book [30] as well as from ([2~42|) . 

Lemma 2.2. Let coefficients of the operator Lq in satisfy conditions il.5]) . 1 1.6]) . Then there 

exist sufficiently large constants vq = Uq (Lq, -RG3/4) > 1, Ao = Ao \Lq, -RG3/4) > 1 such that the following 
Carleman estimate holds 



CAV«p(».4-) / („= + |V„|* + «?)^ 



+CAV H) exp 



I r 



^u 2 + I Vu| 2 + u 2 ^ (icr + / (ut — Lqu) 2 £ 2 {x, t) dxdt 



d 2 D 3 



D, 



>C I (\v\\7u\ 2 + \ 3 v^- 2 "u 2 ^ 2 (x,t)dxdt,Vv>v y\>\ yueC 2A (D 3/i ) 



D3/4 



Theorem 2.2. Let conditions U.3)) - n~B)) and \2.2J$ be valid. Suppose that the function w £ C 2,1 (Q) 
satisfies the following integral inequality 



J (w t - Lw) 2 dxdt < K 2 ,K = const. > 0. 



(2.45) 



Let 



Denote 



(3 (x, t) = w Pi (x, t) = d Xl w g (x) = w (x, 0) , x G O. 



F = \\Po\\m( dl *) + \\Pi\\l 2 ( 9i *) + K - 
9 



(2.46) 



Assume that an upper bound C'2 = const. > for the norm HHIc^q) * s known, 

IIHI C i(*) < C2. 



(2.47) 



Then there exists a sufficiently small number So = So (L, -RG3/4) € (0, 1) dependent only on listed parameters, 
such that if F € (0, So), then the following logarithmic stability estimate is valid 



NhC 2 



(2.48) 



In particular, in the case of IP2, let Assumption holds and \1. 7| ) be valid. Assume that, in addition to the 
above, A2.12\) is valid, and the number 8 in H2.12]) is so small that C8 € (0, So), where C = C (P, L, <!>, d) > 
is the number from \2. 36\) . Also, assume that for a certain a € (0, 1) the upper bound C3 of the norm 
\\f\\c*+°(m ^ given, i.e. \\f\\ C 2+ a (p\ < C 3 . Then 



l/l 



L 2 (n) 



< 



'In 



cs 



(2.49) 



In addition, there exists a number p — p (Lq, RG31A G (0, 1/2) such that if F € (0, So) , then the following 
Holder stability estimate is valid 



< M 1 C 2 F f '. 



(2.50) 



Proof. In this proof £0 = £0 [L, RG3/4) € (0, 1) denotes different sufficiently small numbers associated 
with Lemma 3.1. By (|2~24|) . (pT28]l , (|23T| and (f2733|) 



n c i?G 1/2 c 



Using (225), (|2~^5j) and (ggS]) , we obtain 



y (w t - Lu>) 2 



dxrfi < -ftT exp 



2-4" 



,V£G (0,£ ) 



(2.51) 



(2.52) 



»3/4 



On the other hand, using Lemma 3.1, we obtain for all v > vq,£ € (0,£q) 



(wt — Lw) 2 tp 2 dxdt > I {wt — Low)* ip z dxdt — M\ I ( \w\* + w £ ) ip z dxdt 



2 .„2, 



2 \ ,„2, 



G 



3/4 



'3/4 



»3/4 



> 



exp 



Fix a number v > v$ (Lq, RG3/4) > 1 such that 



C J |Vw| 2 + ^^"V) if 2 (x, t) dxdt - Mj y (ju;| 2 + w 2 ) (p 2 dxdt 

G3/4 

(ll/?oll h 1 ^!*) + ll/^i IL 2 (di*)) (2.53) 

§(!)] y ( w 2 +ivh 2 +- 2 )^. 



G3/4 

c^ 3 



— exp 

£ V £ 



2-4" 



CV 3 /4 



92G3/4 



5V 1 
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(2.54) 



If necessary, decrease e = e (L,RG 3 /4) £ (0, 1), so that Mi < Gvj (2s) ,Ve G (0, £o) ■ Then (|2.53|) leads to 
the following estimate for all e G (0, eq) 

C 



(u>t — Lit;) 2 (p 2 dxdt > — / ^| Vui| 2 + w; 2 ^ <p 2 (x, t) e?ire?i 



'3/4 



G3/4 

c 

£ 

c 



2 ■ 4 l 

£ 

2 /4 



(||A>I 



ll/?i" 2 



£2(91*) 



£ V3 



(w 2 + |Vw| 2 + iv^j da. 



d 2 G 3 /4, 



Combining this with (|2.46p . ()2.47[) and (|2.52l) and decreasing Eq, if necessary, we obtain 



\Vw\ 2 + w 2 j tp 2 (x,t)dxdt < 

G3/4 

On the other hand, by (|2~29|) 



2 • 5" 



F 2 + CC\ exp 



2 /5 
£ 1 3 



,Vee(0,e ). (2.55) 



(\Vw\ 2 + w 2 ) (p 2 (x,t)dxdt > I (\Vw\ 2 + w 2 ^)tp 2 (x,t) dxdt > exp <2 —— I (|Vw| 2 + w 2 ^ dxdt, V£ £ (0, e ) 



'1/2 



Combining this with (|2.55[) . we obtain 



j (|Vw| 2 + w 2 ) tfecft < Cexp (^y-) F 2 + CCl exp 

G V 2 



,V£G (0,£ ) 



Hence, using (|2.54p . we obtain 

J (\Vw\ 2 + w 2 ) dxdt < Cexp 



2 • 5" 



"T ) 



(2.56) 



Gi 



By (j2T2"5]) and (gUgJl G 1/2 C {i G (e/2, 3e/2)} . Hence, the mean value theorem, (f2"3T1) and ([231)1) imply that 
there exists a number t* £ (e/2, 3e/2) such that for all e £ (0, £q) 



1 



II w ( 2; ^*)IIl2(o) < - II w (- t ^*)IIhi(hg 1/2 ) ^ Cex P 



2 • h v 



F 2 + CC 2 exp ( - — 



(2.57) 



We have 



it) (x, t*) = g (x) + J wt (x,T)dr. 


Hence, using (|2.47p . we obtain 

\\w{xX)\\l 2( n) > Hflll£ a (n)-e|KllL(*) ^ NlL(n) ~ ^f£. 
Combining this with (|2.56p . we obtain 

||.9lli 2( o) < MiC 2 e + Mi exp ^ 2 + MiCg exp (-^J , Ve G (0, e ) • 
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(2.58) 



Choose e — e (F) such that 



cxp ( — — ] F 2 = cxp ( -— I . (2.59) 



e 



Hence, 

e =^±^ ia = 2.&+r. (2.60) 

To ensure that e is sufficiently small, i.e. e £ (0,£o) , we need to choose F so small that 

2 



< F < exp - 

V £ o 

Hence, we choose 5q — 5 (L,i?G 3 / 4 ) = exp (— 2/e ) ■ Hence, (|2.58j) . (|2.59|) and lead to 

MU s + * (' + <*) (O 1 " - 5^J + M ' <' + <*> ( f """) "° • (2 ' 61 » 

as long as A S (0, <$o) ■ Decreasing, if necessary Jo, we obtain (|2.48|) from (|2.61[) . 

We now prove (H3SJ). Let the function v £ c 2+a - 1+a ' 2 (R" x [0,T]) ,VT > be the solution of the 
problem (23), (|23]l . Recall that by dUTJ) / (x) = in R n \ft. It follows from the formula (14.6) of §14 of 
Chapter 4 of the book [21] as well as from (|2.33[) that 

IMI c i(5) < ll«|| oa+ -.i+«,/>(5) < M i Wf\\c^(n) < ( 2 - 62 ) 
By ([2~4] A" = in (|!£4"5|) . Next, since (j2~31))) follows from (|2"7T2]) . we obtain for the new number F in (|2T35|) 

^~ll^ll^(^) + p2|L 2( a 1$) <^- ( 2 - 63 ) 

Thus, (l2~47l) and ([2~48)l imply that (j2~49|) follows from (f2T62|) and ([2~63]) . 

Finally, we prove (|2.50p . Using (|2.37p ~ (|2.44[) as well as Lemma 2.2, we obtain similarly with (|2.50l) 

J (|Vw| 2 + uj 2 ) dxdt < Cexp(2A ■ b v ) F 2 + CCfexp (-2 V+1 A) ,VA > Ai, 

where Ai = Ai (L, RG3/4) > 1 is a sufficiently large number. Hence, we obtain similarly with (|2.59|) and 
(I2~6TD 



J (\Vw\ 2 + w 2 ) dxdt < M x ClF 2p . □ 



3. The Quasi- Reversibility Method (QRM). We construct the QRM only for the more difficult 
case of IP 2. The case of IP1 is similar, and it can be derived from [TB]. Also, we work in this section only in 
3-d, keeping the same notations as above. The construction in the 2-d case is similar. Since it was described 
in subsection 2.1 how to stably obtain the Neumann boundary condition in the parabolic case for both IP1 
and IP2, we assume now that we have both Dirichlet and Neumann boundary conditions at 

v |&*= <^ 2 fa, *) ! 9 Xl v | 9l $= V 2 fa, t) . (3.1) 
The QRM means in our case the minimization of the following Tikhonov functional 

J 1 (v) = \\v t - Lv\\ 2 L 2($) + 7 ||u||ff4 (#) , (3.2) 
12 



subject to the boundary conditions (|3.ip . In (|3.2[) 7 > is the regularization parameter, which should be 
chosen in accordance with the level of the error in the data. 

The requirement v G H 4 ($) is an over-smoothness. This condition is imposed to ensure that v G C 1 (W\ : 
because of (|2.47j) and the embedding theorem. However, the author's numerical experience with the QRM has 
consistently demonstrated that one can significantly relax the required smoothness in practical computation, 
see [TH [22] [23] and chapter 6 of [5] . This is likely because one is not using an overly small grid step size 
in finite differences when computing via the QRM. Hence, one effectively works with a finite dimensional 
space with not too many dimensions. This means that one can rely in this case on the equivalence of all 
norms in finite dimensional spaces. Thus, most likely one can replace in real computations 7 |MI#4($) with 

tI|w|Ih2,i ($) . 

While (|3.2j) is good for computations, to prove convergence of the QRM, we need to have zero boundary 
conditions at <9i<£>. Assume that both functions Tp 2 , tp 2 G H 21 (d\<&) . Denote 

r (x, t) = Tp 2 (X, t) + X\l\) 2 (x, t) = Tp 2 (x, t) + X!lp 2 (x, t) , 

v (x, t) = v (x, t) — r (x, t) ,p(x,t) = — (r t — Lr) [x, i) , 
f(x)=v(x,0)=f(x)-r(x,0). 
Using p3J), (|2~5l) and (|3~Tj) . we obtain 

Vt — Lv = p (x, t) , {x, t) G (3.3) 
v I o 1 <s>^0,v Xl \d^= 0. (3.4) 

Thus, we have obtained Inverse Problem 3. 

Inverse Problem 3 (IP3). Find the function f (x) for x £ fl from conditions iS.3]) , {3.$ - 
To solve IP3 via the QRM, we minimize the following analog of the functional (|3.2[) 



J 1 (v) = \\v t - Lv-p\\ L2m +j\\v\\ H 4 m ,veH$($), (3.5) 
i? 4 ($) : = {u e H 4 ($) : u \ M = u Xl |a 1$ = 0} . 

Let (,) and [,] be scalar products in L2 ($) and H 4 ($) respectively. Let the function u 7 € Hq ($) be a 
minimizcr of the functional (|3.5I) . Then the variational principle implies that 



(dtUj — Lu, dtw — Lw) + 7 [u, w] = (p, wt — Lw) , £ Hq ($) . (3-6) 

Lemma 3.1 follows immediately from the Riesz theorem and p.6[) . 

Lemma 3.1. For every function p £ Li ($) and every 7 > there exists unique minimizer it 7 = 
u 7 (p) £ Hq ($) of the functional A3.5\) . Furthermore the following estimate holds 

Mi 

IKIIfl*(») < lbll i2 ($) • 

The idea now is that if u~ (x, t) € Hq ( < 3 ) ) is the minimizcr mentioned in Lemma 3.1, then the approximate 
solution of IP3 is 

/ 7 (x) = u 7 (x,0) . (3.7) 

The question of convergence of minimizers of J 7 to the exact solution is more difficult than the existence 
question of Lemma 3.1. To address the question of convergence, we need to introduce the exact solution as 
well as the error in the data, just as this is always done in the regularization theory [31 [SJ [37] ■ We assume 
that there exists an "ideal" noiseless data p* G L 2 (<&). We also assume that there exists the ideal noiseless 
solution v* G Hq (<&) of the following problem 

— Lv* = p* (x, t) , (x, t) G (3.8) 
tJ* I 9l $ = 0,v* Xl | ai$ =0. (3.9) 
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Let uj £ (0, 1) be a small number, which we regard as the level of the error in the data. We assume that 

||p-P*lli a (*)<<"- ( 3 - 10 ) 

Remark 3.1. For brevity, we work in this section with the parabolic IP3. Still, Theorem 3.1 can 
be easily linked with the original hyperbolic IP2. Indeed, to ensure that p £ L 2 ($) ,we need Tp 2 ,ij) 2 £ 
H 21 (<9i$) . While Assumption implies (|2.7[) . which, in turn guarantees that Tp 2 £ H 2 ' 1 (c?i<&) , there is 
no guarantee that ip 2 "= # 2,1 (di®) (see (|2.10[) V To ensure the latter, we should replace in Assumption 
C 3 (R n x [0,T]) with C* 5 (R™ x [0,T])\ In this case (j2~T0| would be replaced with p 2 || c :>+<«.i+«/2(p I ) - 
C (P, L) ||^2llc 4+Q - 2+Q / 2 (Pi) ' l a ttcr means, in turn that the comparison of the functions p with the exact 
function p* in (|3.10[) would be replaced with the comparison of the approximate and exact data ip 2 and ip 2 
of IP2. This can be done via a routine procedure by replacing (|2.12[) with \\<p 2 — ^2 1 1 c* 5 (p t ) — $ ebT > 0. 
In this case we would have in (|3.10[) uj = u> (S) . 

Theorem 3.1 establishes the convergence rate of the QRM. Note that an upper estimate of the exact 
solution is often assumed to be known in the regularization theory, also see Remark 2.2. 

Theorem 3.1. Let conditions \2.21$ , 13.3\) . I{3.4\ ) and i3. 1 0\) be satisfied and the regularization parameter 
7 in S3. 5)) is chosen such that 7 = 7 (uj) = uj £ (0, 1). Let the function it 7 £ Hq ($) be the unique minimizer 
of the functional S3. 5)) . which is guaranteed by Lemma 3.1. Let the upper estimate Y — const. > for 
the exact solution v* £ Hq ($) be known, \\v*\\ Hi ^ < Y. Then there exists a sufficiently small number 
ujq = luq (L, $) £ (0, 1) such that if uj is so small that uj^J (Y 2 + 1) £ (0,ujq) , then the following logarithmic 
convergence rate takes place 



r-f. 



M{Y 

< , (3.11 

L 2 m ~ ^ln (uj- 1 ) 



where the function f 7 t n \ (x) is defined in \3. 7\ ) and f* (x) = v* (x,0) . Ln addition, for every uj £ (0,Wo) 
there exists a number p — p (L, $) £ (0, 1/2) such that the following convergence rate takes place 

V-^W^^^M^". (3.12) 

Proof. It follows from (|3.8p and (|3.9p that the function v* satisfies the following analog of (|3.6p 

(v* t ~ Lv*,w t - Lw)+~/[v*,w] = (p,w t - Lw) +j[v*,w],Vw £ Hfi ($) . (3.13) 
Let v = u 1 — v * £ H% ($) and p = p - p* £ L 2 ($) . Subtracting (|cTT5)) from flM}, we obtain 

(v t — Lv, u>t — Lw) + 7 [v, w] ~ (p, wt — Lw) — 7 [v* , w] , Vu> £ Hq (<I>) . 
Setting here w := v and using Cauchy-Schwarz inequality and (|3.10p , we obtain 

{v t -Lvfdxdt + 1 \\v\\ 2 Him <uj 2 + 1 \\v*\\ 2 Rim < uj 2 + 1 Y 2 . (3.14) 

Since 7 (uj) = uj £ (0, 1) , then (|3.14[) implies that ||u|| ff 4($) < Y + 1. Hence, using again (|3.14p as well as 
embedding theorem, we obtain with the constant c = c (<I>) > depending only on the domain $ 

INI C i(5) <<*, (3-15) 
(v t - Lvf dxdt < (Y 2 + 1)uj 2 , (3.16) 



/ 



We now apply Theorem 2.2. Comparing (|3.16p and p. 151) with (|2.45p and (|2.47p respectively, we set 



K := F := uj^(Y 2 + 1), C 2 := cY. (3.17) 
Therefore, (j3~TTj) and (|3~T2l) follow from (|3~T7)) . (|2"^g)) and (f^3U]) . □ 
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